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Abstract 


We present a macroscopic anatomical behaviour of skeletal muscles. According to the data available in the literature, we propose 
a first model for the muscle isometric force generation depending on the muscle length and neuronal excitation frequency. We 
moreover provide a more physical-inspired model of the isometric force exploiting a nonlinear spring description with controllable 
characteristics, such as stiffness and rest length. Then, a hybrid model is proposed to describe the sensory system that provides the 
Central Nervous System (CNS) with the measurements of muscle length and its rate of change. 
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1. Introduction 


Skeletal muscles, or simply muscles, are organs of the mus- 
cular system whose cells have the characteristic to produce force 
and movements. They are mostly attached to the bones of the 
skeleton by tendons, and thanks to their contractile and elas- 
tic properties they allow the relative motion of the bones of the 
involved articulations. Thus, the controlled contraction of the 
muscles allows for the performance of everyday movements. 

Each skeletal muscle consists of thousands of muscle fibers 
wrapped together by connective tissue sheaths, see Fig{i] The 
individual bundles of muscle fibers in a skeletal muscle are 
known as fasciculi. Each muscle fiber is comprised of several 
myofibrils containing many myofilaments. When bundled to- 
gether, all the myofibrils get arranged in a unique striated pat- 
tern forming sarcomeres which are the fundamental contractile 
unit of a skeletal muscle. The two most significant myofila- 
ments are actin and myosin filaments which are arranged dis- 
tinctively to form various bands on the skeletal muscle [1][Ch. 1], 
and [B]. 

Skeletal muscles show up with neuronal innervations from 
the sensory-motor system. A set of neurons, the @-motoneurons, 
constitutes the neuromuscular junction and are responsible, ac- 
cording to their excitation frequency œ measured in pulses per 
second, for sarcomeres contractions. And hence they allow the 
generation of the muscle active force. On the other hand, an- 
other set of neurons, constituting the proprioceptive receptors, 
carries from the muscle to the CNS the values of muscle length 
and rate of change (/a- and I- Spindles), and the measurement 
of the force value at the tendon level (Jb-Spindle or Golgi Ten- 
don Organs (GTO)). 

We can thus see the muscle complex, from a ‘System En- 
gineering’ standpoint, as a complex Multi-Input Multi-Output 
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system with the actuator (@- motoneurons), the plant (the mus- 
cles cells, whose contraction produces the exerted force on the 
joint), and the sensors (Ja- and JI- Spindles and GTO), and the 
CNS as a feedback controller. In FigRja schematics of the in- 
volved connections in reported. 

The modelling of how muscles exert force started with the 
work of Hill in which the muscle description via a spring- 
like behaviour is presented. In particular, the relationship takes 
into account a contractile element that produces muscles force 
and this is in series to a passive element representing the ten- 
don model between the muscles fibers and the connected bones. 
This model has then been exploited to generate more detailed 
models [6][7] in which they describe the microscopic behaviour 
of contraction behaviour due to the effects of sarcomeres and 
cross-bridge. In [8], Hill’s model has been enlarged by adding 
a ‘damping element’. In [9] [10| [I], experimental data has 
been exploited to construct an analytic expression of the total 
force generated by the muscles. In these works, the total force 
has been split into length-dependent force Fy, and a velocity- 
dependent force Fy, both depending on the &-motoneuron ex- 
citation frequency. The Fz, term is then split, similarly to Hill’s 
model, into an active plus passive term. The active force de- 
scribes the part of fibers contraction when excited and produces 
movements, while the passive force only refers to the elastic 
mechanical behaviour of the muscle fibers. 

Experimental data and models of isometric (fixed muscle 
length) and isokinetic (fixed muscle velocity) forces are avail- 
able and presented in [9] {10} [11] [12] [13]. 

A system engineering perspective of the neuromuscular sys- 
tem has been presented in in which they give a 
control-oriented interpretation of the involved parts. In 
we find a work dealing with the optimal co-activation problem 
while optimizing the mechanical impedance of an antagonist 
actuated biological joint (such as an elbow). On the other hand, 
in the authors aim at answering how much reduced the 
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Figure 1: Muscle structure. 


muscle model must be to be able to describe the neuromuscular 
system. In recent years, with the development of the contrac- 
tion theory in the control field, the contractile properties of the 
neuromuscular system (involving Hill’s model equations) have 
been studied in [18]. In an analytic model of complete 
and incomplete tetanus contraction has been proposed, and so 
far, according to the authors’ knowledge, gives the closest to 
reality, although incomplete, dynamical model of the muscle 
behaviour. 

The sensory part involved in the muscle system is related to 
the measurement of the plant variables, such as force at the ten- 
don level (GTO), and muscle/fibres length and velocity (Ia- and 
II-spindles). In the last century, they have been studied in depth 
and the most important works have been collected in 
BI], see for a more recent review on the topic. From a sys- 
tem engineering perspective, the GTOs dynamical model has 
been proposed in by exploiting a proper transfer function 
between the force exerted on the tendons and the GTO neurons’ 
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Figure 2: Muscle and Central Nervous System schematics of the connections. 
Picture taken from {4}. 


firing frequency. For Spindles Ia and II, in [24], a static model 
has been proposed to describe the relationship among the mus- 
cle/fibres length, velocity, and the y-motoneurons excitation to 
the spindles firing rate. Such y-motoneurons are in charge of 
modulating the firing responses of the connected Spindle sys- 
tem. 

To the best of the authors’ knowledge, there is no dynamical 
model of the Spindles organs available in the literature. 

The sensory part plays an important role in the movement 
control architecture. In particular, the presence of internal loops 
between the muscles and the spinal cord allows the possibility 
to actuate the reflexes [26], such as: stretch reflexes, 
reciprocal inhibition, and others. The spinal cord, in this setup, 
is in between the CNS and @-motoneurons via the interneuron 
and it is unavoidably involved in movement control [27)[28] 29]. 

The paper is organized as follows, after having provided 
motivations for the work, in Section [2] we give some authors’ 
considerations about the available data. We give in section3] an 
overview of the most relevant models available in the literature. 
Then, in section [4] we describe the developed models for the 
isometric muscle force. Eventually, before the conclusions in 
sec|6| we provide in sec|5] a hybrid modelling of the sensory 
system available in the muscles. 


1.1. Motivations 

This work aims to describe neuromuscular behaviour in a 
control theory paradigm. In particular, we would like to better 
understand the purposes of the neuromuscular loops and how 
they are exploited in controlling limbs movements, see for ex- 
ample for the application of a stabilizing PID action on a 
single-link biomechanical model or to see how this avail- 
able measurement feedback can be exploited in position and 


movement control. The reader can furthermore see to have 
an insight on how waveform signals can be used to control 
delayed closed-loop systems, thus mimicking the spiking be- 
haviour of biological nervous systems. 

We strongly believe that this modelling idea would be a 
good starting point to eventually comprehend the purposes of 
the brain parts involved in the learning and control tasks of dex- 
terous movements, see for example to have 
a physiological and a biomechanical point of view of brain parts 
involved in movements control. 

The final objective will then be to have a complete picture 
and understanding of the structure used by the brain in con- 
trolling human dexterous movements and thus mimic such a 
strategy in the field of nonlinear control systems such as non- 
linear output regulation as in [37]. The study interleaves with 
the concept of inverse and forward models introduced, in neu- 
roscience, by Wolpert and Kawato in their works [39], see 
also for another application of these models. So an inter- 
mediate step of the study will be to better understand how these 
models are formed, learned, and updated in the brain and how 
they are involved in the control strategy. 


2. Preliminaries on available experimental data in litera- 
ture 


In the muscle force, we can distinguish among three main 
components: a passive force F, that is always present due to 
the natural behaviour of tissues and only depends on the muscle 
length, we then have the active static (isometric) force Fz only 
depending on the muscle contraction (length L) and an active 
dynamic (eccentric and concentric) force Fy also depending on 
the rate of change of muscle contraction L. The last two compo- 
nents are active terms and indeed they depend on the activation 
frequency & of the @-motoneuron. 

Data of isometric force available in literature are mainly 
those in the figures of [10], [11], and [16], reported in this sec- 
tion for completeness. It is worth to notice that such data only 
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Figure 3: Normalized Isometric Force. Relationship between the static Force 
and the muscle length at different activation frequency. Plot taken from : 
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Figure 4: Isometric Force. Relationship between the static Force and the mus- 
cle length at different activation frequency. This plot is half useful due to the 
unknown length, but we can notice that there are different rest length (length at 
which the force is null) depending on the activation frequency. Plot taken from 
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Figure 5: Normalized Isometric Force model with correction term given by a 
negative passive force in . Here the activation frequency is normalized with 
respect to its maximum value @ = 120pps 


per second), because 120pps is the maximum excitation fre- 
quency the @-motoneuron can ‘apply’ on the muscle fibre. 

At a = 120pps the active force Fy (i.e. the total force with- 
out passive F, and dynamics Fy terms) has a maximum value 
indexed in literature by the symbol Fj. This maximum value 
F, is associated to a specific length value L, (which is not the 
maximum length value, but the length L corresponding to the 
maximum force value F, at the maximum excitation frequency 
a = 120pps). Almost all available data and proposed models 
in the literature are normalized with respect to these two val- 
ues. Indeed we find force and length data of the form F /F, and 
L/L, respectively. Due to this fact, we also consider a normal- 
ized model for such normalized values. 

A crucial point in our model development is the presence in 
the force characteristic of rest length values (i.e. values of L for 
which the exerted force is null) that depends on excitation rates 
a, see Fig/4] And this fact leads us to consider a rest length 
function Lo (which is not L, as above) depending only on the 
excitation frequency œ and this might be used to describe the 
contracting behaviour of the muscle. A major problem is that 
most of the data available in the literature lack these rest-length 
data which is fundamental for our model development. We thus 
exploit the same technique introduced in [10]. In practice, they 
first consider a model without the presence of rest length (based 


on Gaussian terms, as it will be clearer at the end of section 3}, 
and then they propose to add a negative passive term to compen- 
sate for the lack of rest lengths values in the model and obtain 
the rest lengths points saturating the force at zero. The approach 
is very clear in Fig. 

In the next section, we provide an overview of the two mod- 
els available in the literature. 


3. Literature models comparison 


In our opinion, the most relevant and complete models in 
the literature are those available in |11] and [10]. Here we give 
a short description of such models for the sake of completeness. 


3.1. Brown’s et al model 


In [1], as the Hill’s model [5], the authors consider both 
the active isometric and the dynamic forces generators gath- 
ered inside a contractile element Fcg along which passive term 
Fpg describes the passive force (previously indicated as F,,, but 
we decided to stay faithful to the original nomenclature). The 
resulting total force Fror presents as 


Fror = Fee + Fee (1) 
and For is expressed as 
Fog (o,L,L) = R(a)A(a,L)Fz(L)Fv(L,L) (2) 


where R is a fibre recruitment percentage and describes the 
amount of fibre involved in the force generation depending on 
the activation frequency, A is the term relating the muscle acti- 
vation to the excitation frequency œ and it would have values 
between 0 and 1 (the same values range has been considered 
for R), F; and Fy describe respectively the force-length rela- 
tionship and the force-velocity, and their product provides the 
force expression of the muscle depending on the current length 
L, the velocity of contraction L, and excitation frequency æ. 

In their work, a simplified version of the force function is 
provided fixing R(@) = 1 and making assumptions of quasi- 
statics (close to isometric) conditions during motion and of a 
starting length equal to the rest length, they hence move to a 
force expression of the type: 


Fola, L,L) =Ar(a,L) (FL(L)FV(L.L) +F,(L)) (3) 


where 
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In this model, one can notice that the length Fz and velocity 
Fy dependent force terms are multiplied one each other. This 
might lead to the interpretation that the velocity terms Fy are 
only acting as scaling functions on the ‘supporting’ isometric 
force Fz. We rather prefer to split the isometric Fz and dy- 
namics Fy terms to analyse their contribution individually in a 
spring-damper fashion. Moreover, notice that when È — 0* in 
the equation we have Fy # 1, thus the model loses the conti- 
nuity with the static isometric properties of the muscle after a 
lengthening phase. 


3.2. Winters’ Model 


In [10], Winters work goes deep inside the basic muscle el- 
ements. Indeed in its work, we can find a microscopic descrip- 
tion of the muscles components such as intrafusal and extra- 
fusal elements with their relative damping and stiffness terms, 
and these are connected in parallel and series fashion. Of their 
work, we only focus on the isometric force that has been mod- 
elled as 


Fiso(@,Lcr) = Fr (0,Lcr) + Fp(Lce) 


where, for the isometric force model, they exploit a Gaussian 
function with variable mean value 


LCE 2 
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F,(a,Lcp) = ae ( (5) 
where Lmo, in Winter’s notation, is the optimal muscle length 
value Lo (in this work notation), Log = Lm — Lio — xsg (Lmo + 
Lo) and xsg is the (dimensionless) extension of the extrafusal 
series of spring-element, Lmt is the overall musculotendinous 
length and L;o is the tendon rest length. The value of œ is nor- 
malized with respect to its maximum value. The passive force 
is instead given by the expression 
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Notice that, as touched above, the passive force has a negative 
term for TE <1. 

We strongly believe that, for our study purposes, the func- 
tions exploited in the models presented in this section are very 
complicated to manipulate and use in control analysis. More- 
over, the parameters meaning is very difficult to interpret from 
a physical point of view. 


4. Proposed Models 


Our main objective is to describe the isometric force (Fz + 
F,) as a nonlinear controlled spring with variable stiffness K 
and rest length Lo 


Fiso(0,L) = K(t,L) (L —Lo(a)) (7) 
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Figure 6: Data comparison between the available data in and the evalua- 
tion/intermediate model. r indicates the residuals vector. 


where & is the &-motoneuron excitation frequency and L is the 
muscle length. In order to obtain such a final model, we first 
designed an intermediate/evaluation model following Winters’ 
considerations. Indeed, we exploited the negative passive term 
introduced in to obtain the rest length values for the differ- 
ent excitation frequencies. This step allows us to then construct 
and fit the parameters of the final/control model (7). 


4.1. Evaluation Model 


In order to construct the intermediate/evaluation model we 
consider the isometric force as divided in active Fz and passive 
F, terms, as exploited in [10]: 


Fiso(a,L) = F,(a,L) +F,(L). 


Inspired, then, by the modelling approach in we describe 
the active force curves as Gaussian functions with mean value 
L and standard deviation o and amplitude A as a function of the 
excitation frequency œ, as shown in (4.1). 
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By running a curve fitting algorithm we are able to fit the 


Table 1: Parameters of Evaluation Model active force. 


Cil Ci2 C73 Cid Cis 
1.134 0.723 0.0504 | 11.3297 | 643.259 
0.574 —0.0211 0.8477 0.0008 n.a. 

0.4518 | —0.01704 | —0.1235 | 0.04898 n.a. 
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=——®—— Passive data 
—@— Evaluation Model 
Control Model 


= 
w 
r 


Passive force 
O 
ine) 


12 1.25 1.3 1.35 
Length 


1.15 


Figure 7: Passive force comparisons among the available data in green, the 
Evaluation model in blue, and the control model in magenta. The residual val- 
ues for the Evaluation model has norm |rgy|2 = 0.0077 and for the Control 
model the norm is |rcm|2 = 0.1706 


active isometric force with a very good order of approxima- 
tion, i.e. the standard euclidean norm of the residual vector is 
|r|2 = 0.051, where r is the stack of the differences between 
the available force data and the related value of the proposed 
model. In Fig|6]we plotted the active force comparison between 
the available data from and the intermediate model estima- 
tion. The identified coefficients are reported in Tab[i] 

Exploiting the data available in and their negative pas- 
sive terms considerations we propose a unified continuous model 
of the passive force term exploiting a sigmoidal and exponential 
function 


— h-fi fol 
Fol) = fit oan + A (8) 
were the coefficient values f;, for i = 0,...,6, can be found 


in Tab[2| The sigmoidal behaviour in the passive term, see 


Table 2: Parameters of Evaluation Model passive force. 


fi h B fa fs 
0.00381 | 2.16582 | 8.16634 | 1.45388 


Fig{T0] does not have an experimental data correspondence but 
provides a plausible description of the physical limit of every 
elastic material, i.e. beyond the ‘boundary’ length value the so- 
licited fibres start to break. Then the total isometric force will 
result from the saturation at 0 of all negative terms, as in 


Fpassive (L) + F,(a,L) 


9 
0, Fiso <0 a 


Fıso(&,L) = [ 
4.2. Control Model 


Due to the high complexity of the intermediate/evaluation 
model (9), and in order to make the force relationship more in- 


tuitive from a physical standpoint, we exploited the intermedi- 
ate model to find the rest length points and then the equivalent 
stiffness values. These data were fundamental to then fit the 
rest length Lo(@) and stiffness K(a,L) functions’ parameters 
in terms of the œ and L, so to construct the final/control model. 
Thus, the isometric force function is given by (7). reported here 
for completeness, 


Fiso(0,L) = K(at,L) (L—Lo(a)) (10) 


where both Lo and K are taken as polynomial functions in @ 
and L 


Lo(a@) = £07? +hat+ey 


(11) 
K(@,L) = kı (L) a? +ko(L)a? +k3(L)o + k4(L) 
where the k;(-), i= 1...4, a polynomial functions in a 
ki(a) = ki af +k? + kizgat* + kia? + kis 07 
(æ) = ki i2 3 4 5 (12) 


+ kip + kiz. 


The relative coefficients values are reported in Tables Blana [4] 
respectively. 


Table 3: Parameters of Control Model rest legth. 
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We also report in Fig[8]a comparison between the stiffness 
surface K(o,L) and the constructed point from the evaluation 
model with respect to the œ and L values. 


Fitting of the stiffness surface 


Figure 8: Fitting of the stiffness K(@,L) surface with respect to œ and L 


The model provides a good approximation of the Evaluation 
Model and, as a consequence, of the experimental data, as it 
is shown in Fig] Moreover, the new model of the total force 
gives a good approximation of the Evaluation one in the range 
(0.4, 1.5]Lo, see Fig[10] This model still provides a good level 
of approximation of the available data. Indeed, the residual vec- 
tor norm is quite small |r|. = 0.18, especially for L < 1.5. Since 
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Figure 9: Data comparison with Spring Model. r indicates the residuals vector. 


higher length values are usually not admissible in standard mo- 
tions due to the mechanical limitation of biological joints, one 
can exploit the proposed control model in scenarios involving 
such standard motions. 
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Figure 10: Comparison between the Evalution and the Control Model for dif- 
ferent excitation frequencies, i.e. œ € 0,30, 60, 120 pps. 


5. Sensors modelling 


In this section, we talk about the dynamic modelling of the 
organs available as sensors. They measure the muscle length 
and velocity, along with the applied force at the tendon level. 
Respectively, these organs are spindle II and Ia, and Ib (also 
called Golgi Tendon Organs). For a model of the last ones we 
refer the reader to the works and [23]. 


Table 4: Parameters of Control Model stiffness. 
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5.1. Spindles Ia and IT 


Ia- and Il- spindles neuronal activities are correlated to the 
length and rate of change of the muscle fibres contractions. 
Along with each of the Spindles, we find active elements, i.e. 
the y-motoneurons whose excitation frequency modifies the spin- 
dles’ characteristics. In particular, such modulation can be in- 
terpreted as an increasing of sensing precision from the Spin- 
dles organs, despite the noise amplification drawback. We can 
generally distinguish between static, ys, and dynamic, Y4, Y- 
motoneurons. In the modelling we only consider the effects of 
the yz-motoneurons on the spindles spiking frequency due to 
the lack of available data. 

It is worth noticing in Fig[ITb|that in Spindle-Ia during neg- 
ative velocities transient the output of the sensor saturates at 
zero while it is not clear what happens for Spindle-II due to 
the end of y,-excitation during the same part of the trajectory 
(negative velocities transients). 

To construct the model of the sensors, we first find the static 
behaviour of the sensors’ outputs, we make an analysis on the 
slopes as similarly done in [24]. 

We now first show the hybrid dynamical model found for 
the Spindle-Ia and then the one of the Spindle-II. 


5.1.1. Spindle-Ia 
We considered an Hybrid dynamical model and is given by 


Xta = Pla (Kt L — Ris) 

Pla =0 
Kui, = 0 

Yla = Sato (xra + Ky, Yia + Kib) 
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where x7, is the filter state, yj, is its output, the function sato 
is the function that saturates at zero negative values, pja is the 
filter pole, Kz, Ky,, and Ky, are the static gains for the muscle 
length, velocity and the excitation state of the yj-motoneuron; 
in this case, L, V and Yra are the muscle length, velocity and the 
excitation state of y;-motoneuron and are considered as system 
inputs. Then X7, is the vector collecting the state, inputs and 
output of the system. The flow set Cra = R3\Dia while the 
jump set is Dja = Dia U Diaz U Dia3, where 


Drat = {Xma : L < 0& (Kz, #9|| Pia # Pra)}; 


Dia = {Xm : L > O& (Ki, # Ki, \|Pta £ Pta) }, 


64.1912 | —405.645 


120.853 


Yia [PPS] 


3 


mb 


Time [s] 
(a) Spindle Ia simulation with hybrid dynamical system 


oe 
+ 
P 


s 


ferent frequency 


la af 


2.5 


a 5.0 7.5 sec 
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Dia = {Xia : L = O& (Kz, # O|] Pra + Piq)}- 


In Fig{IT]we compare the output of the dynamical system 
with the available experimental data. The identified parameters 
are all gathered in Tab/5] 


Table 5: Parameters of spindle Ia model. 
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(b) Experimental Data of Spindle II in (21)[Ch.19] 


Figure 12: Spindle II Model output with experimental data in [Ch.19] com- 
parison 


5.1.2. Spindle-IT 


Also for this sensing system a Hybrid dynamical model and 
is given, similarly, by 


Xi = pr (oy high = x11) 


Pu =0 
Ci : Cli =0 f 
yu = Sato [x7 + Ky; Yir + Kj,,L+ 
(1~€2,,) Kr, (L—L)] (14) 


(YXr € Din) ői, =9, pi = Pu 
(VX € Din) $k, =l, Pi = Pu 
( 
( 


VX € D3) 6i, =9, Pir = Py 
VX € Dr) xf =X, yf = Y1 

where x,; is the filter state, yyy is system output, the function 
sato is the function that saturates at zero any negative argu- 
ment, pz; is the filter pole, Kz,,, Ky, and Ky, are the static 
gains for the muscle length, velocity and the excitation state 
of yı-motoneuron; in this case, L, V and Yr are the muscle 
length, velocity and the excitation state of yj-motoneuron and 
are considered as system inputs. Then X7; is the vector col- 
lecting the state, inputs and output of the system. The flow set 
Cy = R?\Dy while the jump set is Dy = Din U Din U Din, 


where 
Din = {Xu : L <0& (éL; #0||pu + Pr)}, 
Din = {Xn : L > 0& (r; + 1\|pu + Pr}, 


Dip = {Xn : L = 0& (EL, ZO||pr # Pj,)}- 


In Fig[I2]we make a comparison between the output of the dy- 
namical system with the experimental available. The iden- 
tified parameters are all gathered in Tab|5] 


Table 6: Parameters of spindle II model. 


6. Conclusions and future works 


In this work, we dealt with the modelling of muscle force 
considering the passive and active components of the isometric 
force. We first provided a brief overview of the model available 
in the literature and then show two models, an Evaluation and a 
Control one. The first one is more complex and precise than the 
second but does not add more concerning the available model 
in the literature. On the other hand, it allowed us to construct a 
set of data that have been useful for the identification of the sec- 
ond model parameters. The latter (the Control model) becomes 
easier to manage, although less precise than the intermediate 
one, and gives a physical intuition about the muscle behaviour. 
Indeed, it describes such muscular behaviour as a controllable 
non-linear spring with controllable (variable) stiffness K and 
rest length Lo. 

We furthermore provide a Hybrid dynamical model formu- 
lation of the spindles organs that act as muscle length and ve- 
locity sensors for the Central Nervous System. 

Future works concern the development of the muscle char- 
acteristic depending on the velocity of shortening and lengthen- 
ing of the muscle fibres, so to provide a more complete model, 
and then analyse this new model in a control scenario involving 
joints with antagonist actuation. 
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